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Abstract. The linear inverse source and scattering problems are studied from the perspec- 
tive of compressed sensing, in particular the idea that sufficient incoherence and sparsity 
guarantee uniqueness of the solution. By introducing the sensor as well as target ensembles, 
the maximum number of recoverable targets (MNRT) is proved to be at least proportional to 
the number of measurement data modulo a log-square factor with overwhelming probability. 

Important contributions include the discoveries of the threshold aperture, consistent with 
the classical Rayleigh criterion, and the decoherence effect induced by random antenna 
locations. 

The prediction of theorems are confirmed by numerical simulations. 



1. Introduction 

We consider the imaging problem in the form of inverse source or scattering problem which 
has wide-range applications such as radar, sonar and computed tomography. The imaging 
problem is typically plagued by nonuniqueness and instability and hence mathematically 
challenging. Traditional methods such as matched field processing j25] are limited in the 
number of targets that can be reliably recovered at high resolution. They often fail to detect 
a substantial number of targets, while at the same time they tend to produce artifacts 
obscuring the real target images. These limitations are due to the presence of noise and the 
fact that the imaging problem is in practice under determined. The standard regularization 
methods can handle to some extent the problem with noise but are inadequate to remedy 
the issue of nonuniqueness of the solution. 

In this paper we utilize the fact that in many imaging applications the targets are sparse 
in the sense that they typically occupy a small fraction of the overall region of interest (the 
target domain). This sparsity assumption suggests to approach the imaging problem by 
using the framework of compressed sensing. 

At the core of compressed sensing lies the following problem (here we focus, as is common 
in the compressed sensing community, on the discrete setting). Assume X G is a signal 
that is sparse, i.e., the number of its non-zero components (measured by the ^o-Quasinorm 
||X||o which is simply the number of non-zero entries of X) satisfies s := ||X||o <C m. 
Let y G be the measurement data vector. We explore in this paper the linear inverse 
problem which can be formulated as y = AX where A is an n x m matrix with n <^ m. 
The goal is to recover X, given the data vector Y and the sensing matrix A of full rank. As 
n ^ AX = y is severely underdetermined and unique reconstruction of X is in general 
impossible. 
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However, due to the sparsity of X one can compute X by solving the optimization problem 
(LO) min||X||o s.t. AX = y. 

Since ( [L0| ) is NP-hard and thus computationally infeasible, we consider instead its convex 
relaxation, also known as Basis Pursuit (BP), 

(LI) min||X||i s.t. AX = y 

which can be solved by linear and quadratic programming techniques. The amazing discovery 
due to David Donoho was that under certain conditions on the matrix A and the sparsity 



of X, both (LI) and (LO) have the same unique solution [T^. One such condition is the 
Restricted Isometry Property (RIP) due to Candes and Tao [7J, which requires essentially that 
any nx s submatrix of A is an approximate isometry. This property is satisfied by a number 
of matrices such as Gaussian random matrices or random partial Fourier matrices O EJ [24] . 
In that case, as long as 5 < (9(n/log(m)), with high probability the solution of ( [lT| ) will 
indeed coincide with the solution of ( [L0| ). Another conditon for which equivalence between 
(LO) and (LI) can be proven is based on the incoherence of the columns of A, which refers to 
the property that the inner product of any two columns of A is small p!2l [T71 [26| . Moreover, 
the performance of BP is stable w.r.t. the presence of noise and error [6l[13l[27j. Finally the 
computational complexity of BP can be significantly reduced by using the various greedy 
algorithms in place of the linear programming technique EU [22l [26l EZj . The most basic 
greedy algorithm relevant here is Orthogonal Matching Pursuit (OMP) which has been 
thoroughly analyzed in [26j. 

For the imaging problem, the sensing matrix A represents a physical process (typically 
wave propagation) and thus its entries cannot be arbitrarily chosen at our convenience. 
Therefore we cannot simply assume that A satisfies any of the conditions that make com- 
pressed sensing work. The few physical parameters that we have control over are the wave- 
length A of the probe wave, the locations and number n of sensors and the aperture A of 
the probe array. This is one of the reasons that make the practical realization of compressed 
sensing a challenging task. 

The paper is organized as follows. In Section [2] we describe the physical setup, formulate 
the imaging problem in the framework of compressed sensing and make qualitative state- 
ments of our main results. In Section [3] we prove the main result for the inverse source 
problem, in particular the coherence estimate (Section |3.lD and the spectral norm bound 
(Section 3.2). In Section [4| we prove the main result for the inverse Born scattering problem 



for the response matrix imaging (Section 4.1) and the synthetic aperture imaging (Section 



4.21 ). In Section [5| and Appendix B, we discuss the numerical method and present simulation 
results that confirm qualitatively the predictions of our theorems. In Appendix |A] we discuss 
the RIP approach to our problems. 



2. Problem formulations and main results 

In this paper, we study the inverse source and scattering problems both in the linear 
regime to suit the current framework of compressed sensing. For simplicity and definiteness 
we consider the three dimensional space and assume that all targets are in the transverse 
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plane {z = Zq} and all sensors are in another transverse plane {z = 0}. The exact Green 
function for the Helmholtz equation which governs the monochromatic wave propagation is 



gicj|r-a| 

(1) G(r,a) = — r, r=(x,2/,zo), a=(C,7y,0). 

47r|r — a| 

We assume that the phase speed c = 1 so that the frequency uu equals the wavenumber. 

We consider the Fresnel diffraction regime where the distance Zq between the targets and 
the sensors is much larger than the wavelength of the probe wave and the linear dimensions 
of the domains ^ 

(2) zo^ A + L, zo > A 

where L is the linear dimension of the target domain. This is the remote sensing regime. 
Under ^ the Green function ([T]) can be approximated by the universal parabolic form [3] 

(3) G(r, a) = e^-l--eiV(2.o)g.c.|y-r;|V(2.o)^ 

4:7rzo 

which is called the paraxial Green function. This follows from the truncated Taylor expansion 
of the function |r — a| 

I I ^ \x — ^\'^ ^ \y — Typ 
2zo 2zo 

under g. 

In the case of the inverse source problem, the corresponding sensing matrix A is essen- 
tially made of the paraxial Green function for various points in the sensor array and the 
target domain. In this set-up, the entries ^ of the paraxial sensing matrix have the same 
magnitude so without loss of generality the column vectors of A are assumed to have unit 
^^-norm. 

A key idea in our construction of a suitable sensing matrix is to randomize the locations 
Rj = {0^^j^r]j)^j = 1, ...,n of the n sensors within a fixed aperture (a square of size A for 
example). Indeed, we assume are independent uniformly distributed in [0,^4]. We 

assume that the antenna elements are independently uniformly distributed in a square array 
[0, A] X [0, A] in the plane {z = 0}. Define the sensor ensemble to be the sample space of n 
i.i.d. uniformly distributed points in [0, A]^. 

We consider the idealized situation where the locations of the targets are a subset of a 
square lattice. More precisely, let be a regular square sub-lattice Al = {r^ : z = 1, m} 
of mesh size £ in the transverse plane {z = Zq}. Hence the total number of grid points m is 
a perfect square. We defer the discussion on extended targets to the concluding section. 

Let S = {vj^ : / = 1, 5} be the set of target locations and aj^J = 1, 5 be the (source 
or scattering) amplitudes of the targets. Set = 0,i {ji, ...^js}- Define the target vector 
A to be A = (cTj) G C^. We consider the target ensemble consisting of target vectors 
with at most s non-zero entries whose phases are independently, uniformly distributed in 
[0, 27r] and whose support indices are independently and randomly selected from the index 
set {1, 2, m}. The number s = || A||o is called the sparsity of the target vector. 
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For source inversion the targets emit the paraxial waves described by ^ which are then 
recorded by the sensors. The measurement vector Y can be written as 

(4) Y = AX 
where the matrix A = [Aij] ^ C^^^ have the entries 

(5) Aij = G(a^,r^-), \/i = 1, ...,n, j = 1, ...,m. 

The first main result proved in this paper can be stated roughly as follows (see Theorem 
[2] and Remark [3] for the precise statement). 

Result A. Suppose 

(6) 1^--^^^- 

For the product ensemble of targets and sensors, sources of sparsity up to (9(n/(lnm)^) can 
be exactly recovered by BP with overwhelming probability. 

When only the sensor ensemble is considered, all sources of sparsity up to 0{^/rl) can be 
exactly recovered by BP and OMP with overwhelming probability. 

The relation ^ indicates the existence of the threshold, optimal aperture given by Xzq/I 
corresponding to p = 1 (see Remark [2] for more discussion on this point). Since the meshsize 
£ has the meaning of resolution, p = 1 is consistent with the classical Rayleigh criterion [3j 

(7) £>^- 

Our numerical simulations (Figure [T]) indeed indicate that ([T]) is sufficient to realize the 
performance stated in Result A. 

Next we consider two imaging settings where the targets are scatterers instead of sources. 
For point scatterers of amplitudes aj^ located at r^^ / = 1,2,3, ...5, the resulting Green 
function G, including the multiple scattering effect, obeys the Lippmann-Schwinger equation 

s 

G(r,a^) = G(r,a^) + ^a^-^G(r,r^JG(r^-^,a^), i = l,...,n. 

1=1 

The exciting field G(rjp a^) is part of the unknown and can be solved for from the so called 
Foldy-Lax equation (see e.g. [15] for details). 

Hence, the inverse scattering problem is intrinsically nonlinear. However, often linear 
scattering model is a good approximation and widely used in, e.g. radar imaging in the 
regimes of physical optics and geometric optics O [8j (see [29j for a precise formulation of 
the condition). 

One such model is the Born approximation (also known as Rayleigh-Gans scattering in 
optics) in which the unknown exciting field is replaced by the incident field resulting in 

s 

(8) G{r,ai) - G{r,ai) = ^aj^G{r,rj^)G{rj^,ai), i = l,...,n. 

1=1 



The left hand side of ([s]) is precisely the scattered field when the incident field is emitted from 
a point source at a^. The Born approximation linearizes the relation between the scatterers 
and the scattered field. The goal of inverse scattering is to reconstruct the targets given the 
measurements of the scattered field. 

For the response matrix (RM) imaging \T5^, 16], we use the real array aperture as in the 
inverse source problem discussed above except the array is also the source of n probe waves. 
One by one, each antenna of the array emits an impulse and the entire array receives the 
echo. Each transmitter-receiver pair gives rise to a datum and there are altogether data 
forming a datum matrix called the response matrix. These data represent the responses of 
the targets to the interrogating waves. 

From ([s]) we see that the corresponding sensing matrix has the entries 

where / is related to i, A; as 

I = i{n-l) + k. 

In the second setting, called the synthetic aperture (SA) imaging, the real, physical array 
consists of only one antenna. The imaging aperture is synthesized by the antenna taking 
different transmit-receive positions ai,i = 1, ...,n [16]. 

The SA imaging considered here is motivated by synthetic aperture radar (SAR) imaging. 
SAR is a technique where a substantial aperture can be synthesized by moving a transmit- 
receive antenna along a trajectory and repeatedly interrogating a search area by firing re- 
peated pulses from the antenna and measuring the responses. This can greatly leverage a 
limited probe resource and has many applications in remote sensing. The image formation is 
typically obtained via the matched filter technique and analyzed in the Born approximation 

Here we consider a simplified set-up, neglecting the Doppler effect associated with the 
relative motion between the antenna and targets. In this case, the sensing matrix A^^ has 
the entries 

(9) ylg^ = G'(a„r,), z = l,...,n, i = l,...,m. 

In other words, A^j^ = with I = i{n — 1) + i. A crucial observation about SA imaging 
is that 

(10) G2(a„r,;cj) -G(a„r,;2cc;) 

modulo a Zo-dependent factor which does not matter. 

The following is a rough statement for inverse Born scattering (Theorems [7|[8] and Remarks 
[ij [5]) proved in Section |4| 

Result B. (i) For RM imaging, assume the aperture condition 

For the product ensemble of sensor and target, scatterers of sparsity up to (9(n^/(lnm)^) 
can be reconstructed exactly by BP with overwhelming probability. 

When only the sensor ensemble is considered, all scatterers of sparsity up to 0{n) can be 
exactly recovered by BP and OMP with overwhelming probability. 
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(ii) For SA imaging, assume the aperture condition 

(11) 2/peN 

which is weaker than 

For the product ensemble of sensor and target, scatterers of sparsity up to (!)(n/(lnm)^) 
can be reconstructed exactly by BP with overwhelming probability. 

When only the sensor ensemble is considered, all scatterers of sparsity up to 0{^/rl) can 
be exactly recovered by BP and OMP with overwhelming probability. 



As a result of the SA aperture condition (11), the corresponding optimal aperture is 
half of that for the inverse source and RM imaging. In other words, SA can produce the 
qualitatively optimal performance with half of the aperture. This two- fold enhancement of 
resolving power in SA imaging has been previously established for the matched-field imaging 
technique [16] . 

Our numerical simulations (Section [5]) confirm qualitatively the predictions of Result A 
and B, in particular the threshold aperture and the asymptotic number of recoverable targets. 

Currently there are two avenues to compressed sensing [4J : the incoherence approach and 
the RIP (restricted isometry property) approach. When the RIP approach works, the results 
are typically superior in that all targets under a slightly weaker sparsity constraint can be 
uniquely determined by BP without introducing the target ensemble. We demonstrate the 
strength of the RIP approach for our problems in Appendix [A| (see Theorem [TT| and Theorem 



12 for stronger results than Result A and Result B (ii), respectively). However, Result 
B(i) seems unattainable by the RIP approach at present, particularly the quadratic-in-n 
behavior of the sparsity constraint. On the other hand, the incoherence approach gives a 
unified treatment to all three results and therefore is adopted in the main text of the paper. 

3. Source inversion 

Let G(r, a) be the Green function of the time-invariant medium and let G be the Green 
vector 

(12) G(r) = [G(r,ai),G(r,a2),...,G(r,a„)]* 

where t denotes transpose. For the matrix ^ define the coherence of the matrix A by 

^(^^-^#iiG(p.)iiiiG(p,)ir 

The following theorem is a reformulation of results due to Tropp [28] and the foundation 
of the imaging techniques developed in this paper. 

Theorem 1. Let X be drawn from the target ensemble. Assume that 

(in \ — ^ 
81n— j , ee(0, 1) 

and that for p > 1 

(14) ,(2hiY\±iw\i< 1 
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Then X is the unique solution of BP with probability 1 — 2e — s ^ . Here || A||2 denotes the 
spectral norm of A. 

We explain the connection of the theorem with [28j in Appendix [B| 

Theorem 2. Let the target vector be randomly drawn from the target ensemble and the 
antenna array be randomly drawn from the sensor ensemble and suppose 



(15) 
(16) 



— = - e N. 



- 2 



5, > 0. 



then the targets of sparsity up to 
(17) 



s < 



n 



64 In ^ In ^ 

e 



can be recovered exactly by BP with probability greater than or equal to 



(18) 



1 - 2S 



pn{n — 1)^/^ 



1/2 



[1 - 2e - s-P] 



P 



In m — In e 
288^6 Ins' 



Proof. The proof of the theorem hinges on the foUowing two estimates. 



Theorem 3. Assume (15) and 
(19) 



m < 



for some positive S and K. Then the coherence of A satisfies 
(20) /x(A) < A/2X/vn 

with probability greater than (1 — 5^. 

Remark 1. The general lower bound for coherence ^TOi [30] 



m 



n 



n{m — 1) 



< M < 1 



implies that the coherence bound (20) is optimal modulo a constant factor. 

Remark 2. Since the coherence of the sensing matrix should decrease as the aperture in- 



creases and since the analysis in Section \3.1\ shows that the coherence is of the same order 
of magnitude as n~^l^ whenever (15) holds, simple interpolation leads to the conclusion that 
the coherence should be roughly constant for 



(21) 



A> 



corresponding to p < 1. The right hand side of (21), corresponding to p 
optimal aperture. 
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1, defines the 



(23) 



1 - 



Theorem 4. The matrix A has full rank and its spectral norm satisfies the bound 
(22) ||A||2 < 2m/n 

with probability greater than 

pn{n — 1)^/^ \zq 
Remark 3. By the theorems of Donoho, Elad [12j and Tropp [26j, the targets of sparsity 

can be recovered exactly by BP as well as by Orthogonal Matching Pursuit (OMP). 
Theorems\^ and^ imply that with probability greater than 

1-25- ^^"^ ~ ^^^'^ 
of the sensor ensemble, all targets of sparsity 



m 



1/2 



s <-{! + 



can be recovered exactly by BP as well as OMP. 



(24) 



Condition ( |17| ) implies the existence of K such that 

n 



As a consequence (19) and (13) are satisfied with probabihty greater than 1 — 2(^ by Theorem 

El 



(25) 



Now the norm bound (22) implies (14) if 

pins ^ 



21n- 



1 



2s 

^77 - 4eV4^ 



p> 1, 



which in turn follows from (17) and the condition 

, 2m , m 1/1 
^"X^"7^ 96 2^- 



phis 
21n^ 



1/2^ 



Hence for m 3> 5 (and hence n 5) we can choose p m (14) to be 

In m — In 6 

p 



72^6 In 5 ' 

Since Theorems |3] and |4] hold with probability greater than 

pn{n — 1)^/^ 



25 



m 



1/2 



and since the target ensemble is independent of the sensor ensemble we have the bound (18) 
for the probability of exact recovery. 



□ 



3.1. Proof of Theorem [3l coherence estimate. 

Proof. Summing over a^, / = 1, n we obtain 

n 1 ^ 

(26) ^^^A^-Aij = e'^^i^'j^yj-^i-yiyC^^o) _ ^i^iu;{xi-Xj)/zo^ir]iu;{yi-yj)/zo ^ 



Define the random variables X/, Yj, / = 1, n, as 



(27) 
(28) 



Xi = cos [{^i{xi-Xj) + rii{yi-yj))u/zo] 
Yi = sin [{^i{xi - Xj) + r]i{yi - yj))uj/zo\ 



and their respective sums 



Sri 



1=1 



1=1 



Then the absolute value of the right hand side of (26) is bounded by 
1 



(29) 



\Sn + iTn\ < - {\Sn - ^Sn\ + |T„ - ET„| + |E(5„ + iTn)\) . 

n n 



To estimate the right hand side of (29), we recall the Hoeffding inequality 



Proposition!. Let Xi^ ...^ be independent random variables. Assume that Xi G [a^,6/],/ 
1, n almost surely. Then we have 



(30) 



F[\Sn-^Sn\ > nt] < 2exp 



2nH 



2+2 



for all positive values oft. 



We apply the Hoeffding inequality to both Sn and T^. To this end, we have ai = —1,6/ 
1 , V/ and set 



t = Kl^, K >0. 

Then we obtain 

(31) P [n-^ \Sn - IKSnl > K/V^ < 2e-^'/2 

(32) P [n-^ \Tn - ETn\ > X/v^] < 2e-^'/^ 

Note that the quantities depend on Xi — Xj, yi — yj^ i.e. 

Sji{xi Xj, yi Vj)) Tji{xi Xj, yi Vj)- 
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We use (31)- (32) and the union bound to obtain 

(33) P 

L 

< 2(m - l)e-^'/2 

(34) P 

L 

< 2(m - l)e-^'/2 



maxn ^ \Sn{xi - Xj.yi - yj) - ESn{xi - Xj.yi - yj)\ > K/^/n 



maxn ^ |r^(x^ - Xj.yi - yj) - Er^(x^ - Xj.yi - yj)\ > K/^/n 



Hence, if (19) holds for any smaU number 5 > 0, then the right hand side of (33)- (34) is 
less than 6. 

The third term on the right hand side of (29) can be calculated as follows. By the mutual 
independence of and rji we have 



n 



\E{Sn + iT^)\ 



1 

n 

1 

n 



1=1 



1=1 



= |E (e^^^^(^^-^j)Ao^ £ ^^ir]i^(yi-yj)/zo 

since (^/, ry/, / = 1, n are independently identically distributed. 

Simple calculation with the uniform distribution on [0, A] x [0, ^4] yields 

|]g ^^i^lu;{x^-Xj)/zo^ ]g ^^im^{yi-yj)/zo^ I ^ 









1 














sin^ 













(35) 
with 

0^^- = Auj{xi - Xj)/zo, ^pij = Auj{yi - yj)/zo. 
The optimal condition is to choose A such that 

(36) (pij = ^pij G 27rZ, 



under which (35) vanishes. Condition (36) can be fulfilled for an equally spaced grid as is 



assumed here. Let 



£ min \xi — Xj\ min \yi — yj\. 



The smallest £ satisfying condition (|36|) is given by 
(37) 



£=^, X = 2n/u 
A ' 



which can be interpreted as the resolution of the imaging system and is equivalent to the 
classical Rayleigh criterion. 
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In this case, E(S'^ + iT^) = and hence 

/x(A) < V2K/^/n 
with probabihty (1 — 5)^ under the condition (|19|). 



□ 



3.2. Proof of Theorem [4[ spectral norm bound. 

Proof. For the proof, it suffices to show that the matrix A satisfies 

71 

(38) II -AA* 



m 



*-n\\2 



< 1 



where is the n x n identity matrix with the corresponding probabihty bound. By the 
Gershgorin circle theorem, (38) would in turn follow from 

(39) . . . . . 1 



m 



n — 1 



since the diagonal elements of ^AA* are unity. 

Since {^i^r]i)^i = 1, ...,n are uniformly distributed in [0,yl] x [0,yl], ^ ^j^Vi Vj with 
probability one. 

Summing over r^, / = 1, m we obtain 



n 
m 



1=1 



(40) 
Thus, 

(41) 



^iu;(r]i-r]j)(yi+y^i)/zo _ ^iu;(r]i-r]j)yi/ zq 
X 



]_ _ ^iu;{r]i-r]j)i/zo 



n 
m 



1=1 



where we have used the identity 
(42) 

Let 
(43) 



= mm mm 



1 

< — 

m 



sm 



2zo 



sm 



22:0 



e 

sm - 
2 



sm 



/mu;{r]i-r]j)£ 
2zo 



sm 



2zQ 



< 1/2 



which is nonzero with probability one. For i ^ j the random variables 



Xzq Xzq 
11 



have the symmetric triangular distribution supported on [— p~^] with height p = \z^/{lA). 
Note that is an integer by the choice (36). Hence the probabihty that {k > a} for smaU 
a > is larger than 

(1 - 2pa)^(^-^) > 1 - 2pn{n - l)a, P=^ 

where the power n(n — 1) accounts for the number of different pairs of random variables 
involved in (43). 



Using the inequality that 



sinvr/^ >2k, e (0, 1/2), 



(41) and the choice 



1 /n-1 



we deduce with probability larger than 



m 



a 



1 — 2pn{n — l)a = 1 — 



pn{n - 1)^/^ 



m 



1/2 



the decoherence estimate 



p 




implying (39). 



□ 



4. Inverse Born scattering 



In this section, we consider two imaging settings where the targets are scatterers instead 
of sources under the Born approximation ([s]). 

4.1. Response matrix (RM) imaging. For the coherence calculation, we have 

.,g=l 

n 



1=1 



lp=i 



and thus 



In view of (44) and Theorem [3] the following theorem is automatic. 



Theorem 5. Under the assumptions (15) and (19) the coherence of A satisfies 

with probability greater than (1 — 5)^. 

We now proceed to establish the counterpart of Theorem |4} 

12 



Theorem 6. The matrix has full rank and its spectral norm satisfies the bound 
(45) IIA^^II^ < 2m/n2 

with probability greater than or equal to 

pn\n^ - If/^ _Xzo 

Remark 4. ^45 in Remark\^ Theorems^ and^ imply that with probability greater than 

pn^in^ - 1)3/2 



1-25 



of the sensor ensemble, all targets of sparsity 



m 



1/2 



can be recovered exactly by BP as well as OMP. 

Proof. We proceed as in the proof of Theorem |4} As before, we seek to prove 



(46) 



n 



RM* 



m 



< 



n? — 1 



For the RM setting, (26) becomes 

,2 ^ 



m 



X 



^^"^^ ^ 1 — e'^^im^Vk-w-vyy/zo 

where / i{n — 1) + fc, T = i'{n — 1) + k'. 

We apply the same analysis as (26) here. Let 

KVi + Vk -r]k') 



(48) K min min 



^(6 + ^k- & - 



Xzo 



k 



k 



which is nonzero with probability one. For / ^ V the probability density functions (PDF) 
for the random variables 

^(6 + Cfe - 6^ - CfeO KVi + Vk - r]i' - r]k') 

are either the symmetric triangular distribution or its self-convolution supported on [— 2p~^, 2p~^] . 
In either case, their PDFs are bounded by p (indeed, by 2p/3). Hence the probability that 
{n> a] for small a > is larger than 

(1 - 2pa)"'("'-^) > 1 - 2pn\n^ - 1). 
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With the choice 



1 /n2 - 1 

= a 




2 V vn 

we deduce that 

with probabihty larger than 

. o 2/ 2 .^ . pn2(n2- 1)3/2 

□ 

As before, the above estimates yield the following result. 

Theorem 7. Consider the response matrix imaging with the target vector randomly drawn 
from the target ensemble and the antenna array randomly drawn from the sensor ensemble. 



If (fl^y and (16) hold then the targets of sparsity up to 



(49) 



64 In ^ In 

e 



can be recovered exactly by BP with probability greater than or equal to (18). 

4.2. Synthetic aperture (SA) imaging. In view of ([To]), we obtain 

/i(ASAH)=MA(2a;)). 

The following result is an immediate consequence of the correspondence ([9])-(10) between 
SA imaging and inverse source setting. 

Theorem 8. Let the target vector be randomly drawn from the target ensemble and the 
antenna array be randomly drawn from the sensor ensemble. If 

(50) ^ G N 

and ( fJ^ hold then the targets of sparsity up to 
(51) 



64 In ^ In ^ 

e 



can be recovered exactly by BP with probability greater than or equal to (18). 

Remark 5. As in Remark^ conditions ( [7gp and ( f50| j imply that with probability greater 
than 

^1/2 

of the sensor ensemble, all targets of sparsity 



5 < -(1 + 
2^ 



can be recovered exactly by BP as well as OMP. 
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5. Numerical simulations 



In the simulations, we set Zq = 10000 and for the most part A = 0.1 to enforce the second 
condition of the paraxial regime g. The computational domain is [-250, 250] x [-250, 250] 
with mesh-size i = 10. The threshold, optimal aperture according to Theorem [3] is A = 100. 
As a result, the first condition of the paraxial regime ^ is also enforced. Note that the 
Fresnel number for this setting is 

<^ = 360 »l 

indicating that this is not the Fraunhofer diffraction regime and the Fourier approximation 
of the paraxial Green function is not appropriate [3j . 

We use the true Green function ([T]) in the direct simulations and its paraxial approxi- 
mation for inversion. In other words, we allow model mismatch between the propagation 
and inversion steps. The degradation in performance can be seen in the figures but is still 
manageable as the simulations are firmly in the Fresnel diffraction regime. The stability 
of BP with linear model mismatch has been analyzed in [19] for the case when the matrix 
satisfies the Restricted Isometry Property (RIP), see Appendix [A} 

In the left plot of Figure [l} the coherence is calculated with aperture A G [10, 200] and 
n = 100 for the sensing matrices with the exact Green function (red-solid curve) as entries 
and its paraxial approximation (black- asterisk curve). The coherence of the exact sensing 
matrix at the borderline of the paraxial regime with zq = 1000, A = 1 is also calculated 
(blue-dashed curve). All three curves track one another closely and fiatten near and beyond 
A = 100 in agreement with the theory (Theorem [s]), indicating the validity of the optimal 
aperture throughout the paraxial regime. 

Figure [l] (right plot) displays the numerically found maximum number of recoverable 
source points as a function of n with A = 100 by using the exact (red-solid curve) and 
paraxial (black- asterisk curve) sensing matrices. The maximum number of recoverable tar- 
gets (MNRT) is in principle a random variable as our theory is formulated in terms of the 
target and sensor ensembles. To compute MNRT, we start with one target point and apply 
the sensing scheme. If the recovery is (nearly) perfect a new target vector with one additional 
support is randomly drawn and the sensing scheme is rerun. We iterate this process until 
the sensing scheme fails to recover the targets and then we record the target support in the 
previous iterate as MNRT. This is an one-trial test and no averaging is applied. The linear 
profile in the right plot of Figure [l] is consistent with the prediction (17) of Theorem [2| 



To reduce the computational complexity of the compressed sensing step, we use an iterative 
scheme called Subspace Pursuit (SP) [9j. It has been shown to yield the BP solution under 
the RIP P (see Appendix |A| 

In the scattering simulation, we use the Foldy-Lax formulation accounting for all the 
multiple scattering effect ^J. Hence there are two mismatches (the paraxial approximation 
and the Born approximation) in the simulation. 

In the left plot of Figure [2] the compressed sensing image with RM set-up is shown for 
A = 100 and n = 20. The size of the sensing matrix is 400 x 2500 and 35 targets are 
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Figure 1. (Left) The red-solid and black-asterisk curves are, respectively, 
the coherence for the exact and paraxial sensing matrices for zq = 10000, A = 
0.1, n = 100 as a function of aperture. The blue-dashed curve is the coherence 
for the exact sensing matrix for zq = 1000, A = l,n = 100; (Right) the em- 
pirical, maximum number of recoverable sources with |cr| = 1 v.s. the number 
n of antennas for = 100, A = 1 by using the paraxial (black- asterisk) and 
exact (red-solid) sensing matrices. 




Figure 2. (Left) 35 scatterers are perfectly recovered by compressed sensing 
technique with 20 antennas. The red circles represent the true locations of the 
targets. The plot on the right is produced by the conventional matched field 
processing. 



(nearly) exactly recovered. For comparison, the image obtained by the linear processor of 
the traditional matched field processing is shown on the right. In Appendix [C} we outline 
the rudiments of matched field processing. 
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Figure 3. The empirical maximum number of recoverable scatterers (left for 
RM, right for SA) with \a\ = 0.001 v.s. the number n of antennas (or antenna 
locations) for A = 100. The data for n G [10, 30] in the RM plot is fitted with 
the parabola (blue-dashed curve): —16.4950 + 0.1366 * x^. The wavelength is 
0.1 for the RM case. The SA plot depicts the number of recoverable scatterers 
in four settings: paraxial sensing matrix with A = 0.1 (black- asterisk), paraxial 
sensing matrix with A = 0.2 (blue-dashed), exact sensing matrix with A = 0.1 
(red-solid) and exact sensing matrix with A = 0.2 (green-circled) 

In Figure [3] the numerically found maximum number of recoverable scatterers is depicted 
as a function of the number of antennas for ^4 = 100 and for both RM and SA imaging set- 
ups by using the paraxial and exact sensing matrices. Clearly, both curves are qualitatively 



consistent with the predictions (49) and (51). 



6. Conclusions 

In this paper, we have studied the imaging problem from the perspective of compressed 
sensing, in particular the idea that sufficient incoherence and sparsity guarantee uniqueness 
of the solution. Moreover, by adopting the target ensemble following [28j and the sensor 
ensemble, the maximum number of recoverable targets is proved to be at least proportional to 
the number of measurement data modulo a log-square factor with overwhelming probability. 

We have analyzed three imaging settings: the inverse source, the inverse scattering with 
the response matrix and with the synthetic aperture. Important contributions of our analysis 
include the discoveries of the decoherence effect induced by random antenna locations and 
the threshold aperture defined by p = 1 for source and RM imaging and p = 1/2 for SA 
imaging where p = Xzo/{A£). 

In this paper we have considered the localization of point targets and the determination of 
their amplitudes. A natural next step is to consider extended targets. However our approach 
does not extend in a straightforward manner to imaging of extended targets, as can be easily 
seen. Assume that we model an extended target approximately as an ensemble of point 
targets that are spaced very close together. Clearly, this requires the mesh size £ to be so 
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small as to render p 1. To apply our theorems would then require that the aperture 
and the number of antennas increase without bound. Clearly this is not a feasible way to 
image extended targets via compressed sensing. Therefore a somewhat different approach, 
on which we plan to report in our future work, is required for extended targets. 

Appendix A. Restricted isometry property (RIP) 

A fundamental notion in compressed sensing under which BP yields the unique exact 
solution is the restrictive isometry property due to Candes and Tao [7j. Precisely, let the 
sparsity s of the target vector be the number of nonzero components of X and define the 
restricted isometry constant Sg to be the smallest positive number such that the inequality 

{1 - Ss)\\Z\\l < \\AZ\\l < (1 + SMZWl 

holds for all Z G of sparsity at most s. 

Now we state the fundamental result of the RIP approach. 

Theorem 9. [7| Suppose the true target vector X has the sparsity at most s. Suppose the 
restricted isometry constant of A satisfies the inequality 

(52) 5ss + 3(^4. < 2. 

Then X is the unique solution of BP. 

Remark 6. Greedy algorithms have significantly lower computational complexity than linear 
programming and have provable performance under various conditions. For example under 
the condition S^s < 0.06 the Suh space Pursuit (SP) algorithm is guaranteed to exactly recover 
X via a finite number of iterations [9j . 

In this appendix we show that the sensing matrix for source inversion satisfies RIP. This 
can be readily seen by rewriting the paraxial Green function ^ 

for r= (x,2/,Zo),a= (^,^^,0). 

Now the sensing matrix ^ can be written as the product of three matrices 

A = Di$D2 

where 

Di = diag(e^"(^'+^')/('"°)), D2 = diag(e^"("'+^')/(2.o)) 

are unitary and 

Assume without loss of generality that xi = yi = UJ = 0, ...,m — 1 and suppose that 
{ij'}Vj)->3 = 1) •••5^ ^1*6 independently and uniformly distributed in [0,^] x [0,^] with 

cf. 
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The result essential for our problem is due to Rauhut [2 
Theorem 10. [23j // 

^ > Csln^ slum In 



Inn e 



for e G (0, 1) and some absolute constant C , then the restricted isometry condition (52) is 
satisfied with probability at least 1 — e. 

See [5l El] for similar results for sensors located in a particular discrete subset of [0, ^4] x 
[0,A]. 

Since Di and D2 are diagonal and unitary, A satisfies (52) if and only if $ satisfies the 
same condition. 

Theorem 11. Let the sensor array be randomly drawn from the sensor ensemble satisfying 

> Csln^ slum In 



Inn ~ e 
for e G (0, 1) and some absolute constant C ^ then with probability at least 1 — e all source 
amplitudes of sparsity less than s can be uniquely determined from BP. 



From the relationships ([9]), ( |10| it follows immediately that A^^ also satisfies (52) if 

(55) ^ - i 

^^^^ A^o " 2' 



cf. ((SOJ). 

Theorem 12. Let the sensor array he randomly drawn from the sensor ensemble satisfying 
(5^. If 

n 1 
> Cs\v? slum In 



Inn ~ 6 
for e G (0, 1) and some absolute constant C , then with probability at least 1 — e all scatter 
amplitudes of sparsity less than s can be uniquely determined from BP. 

The superiority of the RIP approach, if it works, over that of the incoherence taken in 
the main text of the paper, is that the uniqueness for BP is guaranteed for all targets of 
sparsity at most s and the target ensemble needs not be introduced. Moreover, the stability 
of solution w.r.t. noise is guaranteed [6l[l9]. However, Theorem [7] for the response matrix 
imaging does not seem amenable to the RIP approach. 

Appendix B. Proof of Theorem [T] 

Theorem [l] is an easy consequence from the following two theorems due to Tropp [28] . 

Proposition 2. [28j Let A be a nxm matrix with full rank. Let Ag be a submatrix generated 
by randomly selecting s columns of A. The condition 

(56) 6(pA^2sln(l + s/2))'/' + ;^||A||^<^, p>l 



m 
19 



implies that 

(57) P(||A:A,-I,||2<a)>l-f-^'' 



Proposition 3. [28j Let X be drawn from the target ensemble. If 

(in \ — ^ 
81n— j , eG(0, 1) 

and if the least singular value 

(59) a^in(A,) > 2-1/2^ \S\ = s 

then X is the unique solution of BP (fLJ^^ except with probability 2e. 



First of all, (13) and (14) together imply (56) and (58) with a = 1/2. Moreover, by 
Proposition [2] (59) holds with probability greater than or equal to the right hand side of 
(57). Hence we need only to derive the claimed bound for the probability of the event E 
that X is the unique solution of BP. This follows from the estimate 

F{E) > p(£;|||a:a,-i,||2<2-1)p(||a:a,-i,||2<2-1) 

> (l-2e)(l-(2/5)^) 

> l-2e-(2/s)f. 

Appendix C. Matched field processing 

Matched field processing (MFP) has been used extensively for source localization in un- 
derwater acoustics and is closely related to the matched filter in signal processing. 
The conventional MFP uses the Bartlett processor with the ambiguity surface 

^ G*(r)yy*G(r) 

^<--' ^ l|G(.)||l 

[25] . The Bartlett processor is motivated by the following optimization problem: Maximize 
the quantity 

(61) T^*yy*ty 

subject to the constraint: 

WW = 1. 

The solution 

W = Y/\\Y\\2 

is the weight vector for the matched filter. In the case of one point source of amplitude ai 
located at Xi, 

Y = (TiG(ri) 

hence 

fTiG(ri) 



(62) W 



ki|||G(ri)||2' 
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Extending (62) to an arbitrary field point r by substituting r for ri we obtain the Bartlett 
processor from (61). 

In general, Y is the n-dimensional measurement vector consisting the received signals 
of the array. For inverse scattering in the RM set-up, there are n measurement vectors 
corresponding to n probe signals. The ambiguity surface in this case is the sum of the n 
ambiguity surfaces for the n probe signals. 

In contrast to the conventional matched field processor, the compressed sensing processor 
utilizing the ^^-minimization [HJIH] or various greedy algorithms P, [211 [22l [26] nonlinear. 



References 

[1] A.B. Baggeroer, W.A. Kuperman and P.N. Mikhalevsky, "An overview of matched field methods in 
ocean acoustics", IEEE J. Oceanic Eng.18 (1993), 401-424. 

[2] B. Borden, Radar Imaging of Airborne Targets. Institute of Physics Publishing, Bristol, 1999. 

[3] M. Born and E. Wolf, Principles of Optics^ 7-th edition, Cambridge University Press, 1999. 

[4] A.M. Bruckstein, D.L. Donoho and M. Elad, "From sparse solutions of systems of equations to sparse 
modeling of signals," SI AM Rev. 51 (2009), 34-81. 

[5] E. Candes, J. Romberg and T. Tao, "Robust undertainty principles: Exact signal reconstruction from 
highly incomplete frequency information," IEEE Trans. Inform. Theory 52 (2006), 489-509. 

[6] E.J. Candes, J. Romberg and T. Tao, "Stable signal recovery from incomplete and inaccurate measure- 
ments," Commun. Pure Appl. Math. 59 (2006), 120723. 

[7] E. J. Candes and T. Tao, " Decoding by linear programming," IEEE Trans. Inform. Theory 51 (2005), 
42034215. 

[8] J. C. Curlander and R.N. McDonough, Synthetic Aperture Radar: Systems and Signal Processing, 
Wiley-Intersience, 1991. 

[9] W. Dai and O. Milenkovic, "Subspace pursuit for compressive sensing: closing the gap between perfor- 
mance and complexity," arXiv:0803.0811 , 
[10] P. Delsarte, J. M. Goethals, and J. J. Seidel, Bounds for systems of lines and Jacobi poynomials. Philips 

Res. Repts. 30:3 pp. 91105, 1975, issue in honour of C.J. Bouwkamp. 
[11] D. L. Donoho, "Compressed sensing," IEEE Trans. Inform. Theory 52 (2006) 1289-1306. 
[12] D.L. Donoho and M. Elad, "Optimally sparse representation in general (nonorthogonal) dictionaries via 

minimization," Proc. Nat. Acad. Sci. 100 (2003) 2197-2202. 
[13] D.L. Donoho, M. Elad and V.N. Temlyakov, "Stable recovery of sparse overcomplete representations in 

the presence of noise," IEEE Trans. Inform. Theory 52 (2006) 6-18. 
[14] D.L. Donoho and X. Huo, "Uncertainty principle and ideal atomic decomposition, " IEEE Trans. 

Inform. Theory (2001), 2845-2862. 
[15] A.C. Fannjiang and P. Yan, "Multi-frequency imaging of multiple targets in Rician fading channels: 

stability and resolution," Inverse Problems 23 (2007) 1801-1819. 
[16] A. Fannjiang, K. Solna and P. Yan, "Synthetic aperture imaging of multiple point targets in Rician 

fading media," SI AM J. Imaging Sci.2 (2009), 344-366. 
[17] R. Gribonval and M. Nielsen, "Sparse representation in unions of bases," IEEE Trans. Inform. Theory 

49 (2003), 3320-3325. 

[18] M. Herman and T. Strohmer, " High-resolution radar via compressed sensing," to appear, IEEE Trans. 
Sign. Proc. 

[19] M. Herman and T. Strohmer, "General deviants: an analysis of perturbations in compressed sensing," 
Preprint, Feb. 2009. 

[20] W. Hoeffding, "Probability inequalities for sums of bounded random variables" , J. Amer. Stat. Assoc. 
58 (1963) 1330. 

21 



[21] D. Needell, J. A. Tropp, and R. Vershynin, "Greedy signal recovery review," Proc. 42nd Asilomar 
Conference on Signals, Systems, and Computers, Pacific Grove, CA, Oct. 2008. 

[22] D. Needell and R. Vershynin, "Uniform uncertainty principle and signal recovery via regularized or- 
thogonal matching pursuit". Found. Comput. Math., DOI: 10.1007/sl0208-008-9031-3. 

[23] H. Rauhut, "Stability results for random sampling of sparse trigonometric polynomials," preprint, 2008. 

[24] M. Rudelson and R. Vershynin, "On sparse reconstruction from Fourier and Gaussian measurements," 
Comm. Pure Appl. Math. Ill (2008) 1025-1045. 

[25] A. Tolstoy, Matched Field Processing in Underwater Acoustics, World Scientific, Singapore, 1993. 

[26] J. A. Tropp, "Greed is good: algorithmic results for sparse approximation," IEEE Trans. Inform. Theory 
50 (2004), 2231-2242. 

[27] J. A. Tropp, " Just relax: convex programming methods for identifying sparse signals in noise," IEEE 
Trans. Inform. Theory h2 (2006), 1030-1051. "Corrigendum" IEEE Trans. Inform. Theory (2008). 

[28] J. A. Tropp, "On the conditioning of random subdictionaries," preprint, 2007. 

[29] H.C. van de Hulst, Light Scattering by Small Particles. Dover Publications, New York, 1981. 

[30] L. Welch, Lower bounds on the maximum cross-correlation of signals, IEEE Trans, on Information 
Theory, 20 (1974), pp. 397399. 

Department of Mathematics, University of California, Davis, CA 95616-8633 
E-mail address: f annj iangOmath . ucdavis . edu 

Applied and Computational Mathematics, California Institute of Technology, CA 91125 
E-mail address: yan@acm . caltech . edu 

Department of Mathematics, University of California, Davis, CA 95616-8633 
E-mail address: strohmer@math. ucdavis . edu 



22 



